##### subgroups analyses ######

rm(list = ls())

# Load Packages 
library(tidyverse)
library(haven)
library(labelled)
library(data.table)
library(cregg)
library(ggeffects)
library(survey)
library(gtsummary)
library(performance)
library(car)
# Load Packages 
library(tidyverse)
#library(devtools)
#devtools::install_github("sirpudika/irpanels")
library(irpanels)
library(haven)



# data prep ---------------------------------------------------------------

# Load data
w7_w8_subgroups <- read_csv("/Users/cbrugge/Desktop/CE Paper/w7_w8_subgroups.csv") %>% 
  as.data.frame() %>% 
  rename(obs = `...1`)


# Convert variables to factors
factor_vars <- c("Price", "Functional duration", "Env. production impact",
                 "Energy efficiency", "Recyclability", "Repairability", "product_type",
                 "env_attitudes_dich", "prefs_mandatory_labelling_recy", "prefs_mandatory_labelling_rep",
                 "prefs_mandatory_labelling_duration", "political_ideology",
                 "importance_repairability", "importance_recyclability", "wtp_env_prot",
                 "Income", "Gender", "Age", "ce_category")

w7_w8_subgroups[factor_vars] <- lapply(w7_w8_subgroups[factor_vars], as.factor)


# reorder levels

w7_w8_subgroups$Price <- factor(w7_w8_subgroups$Price, 
                                levels = c("180 / 250 / 470 CHF", "550 / 850 / 890 CHF",
                                           "920 / 1900 / 1650 CHF", "1300 / 3100 / 3200 CHF"))

w7_w8_subgroups$`Functional duration` <- factor(w7_w8_subgroups$`Functional duration`, 
                                levels = c("2 / 2 / 5 years", "4 / 8 / 10 years",
                                           "8 / 15 / 15 years"))

w7_w8_subgroups$`Env. production impact` <- factor(w7_w8_subgroups$`Env. production impact`, 
                                            levels = c("Very low", "Medium",
                                                       "Very high"))

w7_w8_subgroups$`Energy efficiency` <- factor(w7_w8_subgroups$`Energy efficiency`, 
                                                               levels = c("G", "E",
                                                                          "C", "A"))

w7_w8_subgroups$Recyclability <- factor(w7_w8_subgroups$Recyclability, 
                                              levels = c("Not recyclable", "Limited recyclable",
                                                         "Recyclable"))

w7_w8_subgroups$Repairability <- factor(w7_w8_subgroups$Repairability, 
                                        levels = c("Not repairable at all", "Repairable to a limited degree",
                                                   "Repairable"))
# Specify order of levels for ce index
ce_category_levels <- c("Low", "Medium", "High")
w7_w8_subgroups$ce_category <- factor(w7_w8_subgroups$ce_category, levels = ce_category_levels)

# Define choice formula

f1 <- choice ~ Price + `Functional duration` + `Env. production impact` +
  `Energy efficiency` + Recyclability + Repairability

mm <- mm(w7_w8_subgroups, f1, id = ~ PubId, h0 = 0.5, level_order = "ascending")


# ce index & price two-way interaction ------------------------------------



fit <- glm(choice ~  ce_category + Price, family = binomial(), data = w7_w8_subgroups)


ggpredict(fit, terms = c("ce_category", "Price"), vcov_fun ="vcovCL")

theme_set(theme_bw())

mydf <- ggpredict(fit, terms = c("ce_category", "Price"), vcov_fun ="vcovCL")

print(mydf)

plot(mydf, show_ci = TRUE, show_data = TRUE, connect_lines = TRUE) +
  scale_x_discrete(limits = ce_category_levels) +
  labs(x = "Overall product circularity", y = "Probability of product choice") +
  ggtitle("") +  
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "grey") + 
  theme(
    axis.title.x = element_text(size = 15, margin = margin(t = 30)),
    axis.title.y = element_text(size = 15, margin = margin(r = 30)),
    axis.text.y = element_text(size = 15),
    axis.text.x = element_text(size = 15),
    legend.title = element_text(size = 15, margin = margin(b = 10)),
    legend.text = element_text(size=13)
  )





# ce index environmental production impacts two-way interaction -----------


fit <- glm(choice ~  ce_category * `Env. production impact`, family = binomial(), data = w7_w8_subgroups)


mydf <- ggpredict(fit, terms = c("ce_category", "Env. production impact"), ci_level = 0.95, vcov_fun ="vcovCL")

print(mydf)

plot(mydf, show_ci = TRUE, show_data = TRUE, connect_lines = TRUE) +
  labs(
    x = "Overall product circularity", 
    y = "Probability of product choice",
    color = "Environmental impacts \nduring production"
  ) +
  ggtitle("") + 
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "grey") +
  scale_shape_manual(values = 1:3, labels = c("Very low", "Medium", "Very high")) +
theme(
    axis.title.x = element_text(size = 15, margin = margin(t = 30)),
    axis.title.y = element_text(size = 15, margin = margin(r = 30)),
    axis.text.y = element_text(size = 15),
    axis.text.x = element_text(size = 15),
    legend.title = element_text(size = 15, margin = margin(b = 10)),
    legend.text = element_text(size=13)
  )

# ce index & environmental production impacts interaction (levels) --------

w7_w8_subgroups$interaction_env_ce <- interaction(w7_w8_subgroups$`Env. production impact`, w7_w8_subgroups$ce_category, sep = " & ")

mm_interaction_env_ce <- cj(w7_w8_subgroups, choice ~ interaction_env_ce, id = ~ PubId, estimate = "mm")


mm_interaction_env_ce <-mm_interaction_env_ce[order(mm_interaction_env_ce$estimate),]


mm_interaction_env_ce <- mm_interaction_env_ce %>% filter(level == "Very low & Low"| 
                                                            level == "Very low & High" | 
                                                            level == "Very high & Low" | 
                                                            level == "Very high & High")

interaction_labels <- c(
  "Very low & Low" = "Low environmental impact & low circularity",
  "Very low & High" = "Low environmental impact & high circularity",
  "Very high & Low" = "High environmental impact & low circularity",
  "Very high & High" = "High environmental impact & high circularity"
)
p_interaction_env_ce <- plot(mm_interaction_env_ce, vline = 0.5, feature_headers = FALSE,
                             xlim = c(0.15, 0.85)) + 
  labs(x = "Marginal means", y = "Pr(Choice)") +
  geom_point(size = 3.5, position = position_dodge(width = 0.9)) +
  geom_linerange(aes(xmin = lower, xmax = upper), size = 1, position = position_dodge(width = 0.75)) +
  scale_y_discrete(labels = interaction_labels) +  
  ggplot2::theme_grey() +
  theme(axis.title.y = element_blank(),  
        axis.title.x = element_text(size = 16, margin = margin(t = 20)),
        legend.position = "none", 
        plot.title = element_text(hjust = 0.5, size = 16, margin = margin(b = 20)),
        axis.text.y = element_text(size = 16),
        axis.text.x = element_text(size = 12),
        plot.subtitle = element_text(hjust = 0.5, size = 14)) 


print(p_interaction_env_ce)

